The goal of this notebook is to develop a python class for the purpose of counting and labeling flourescent particles captured in images from a prototype for a low-cost medical diagnostic device.
In this notebook we'll consider traditional image processing techniques .. i.e., those that directly operate on an image to extract scientific information without resort to machine learning techniques. We'll save that imoportant discussion for another time.
Computer Vision Textbooks
Nixon, Mark, and Alberto Aguado. Feature extraction and image processing for computer vision. Academic press, 2019. link
Szeliski, Richard. Computer vision: algorithms and applications. Springer Science & Business Media, 2010. [Hesburgh Library][Szeliski Web Page and Materials][Preprint of 2nd Edition]
Programming Books
Howse, Joseph, and Joe Minichino. Learning OpenCV 4 Computer Vision with Python 3: Get to grips with tools, techniques, and algorithms for computer vision and machine learning. Packt Publishing Ltd, 2020.
Villán, Alberto Fernández. Mastering OpenCV 4 with Python: A practical guide covering topics from image processing, augmented reality to deep learning with OpenCV 4 and Python 3.7. Packt Publishing Ltd, 2019.
Pajankar, Ashwin. Raspberry Pi Computer Vision Programming: Design and implement computer vision applications with Raspberry Pi, OpenCV, and Python 3. Packt Publishing Ltd, 2020. Amazon
Papers
Coelho, L.P. 2013. Mahotas: Open source software for scriptable computer vision. Journal of Open Research Software 1(1):e3, DOI: http://dx.doi.org/10.5334/jors.ac
Van Noorden, Richard. "Publishers launch joint effort to tackle altered images in research papers." Nature (2020). https://doi.org/10.1038/d41586-020-01410-9
Blogs and Postings
See https://www.analyticsvidhya.com/blog/2021/04/top-python-libraries-for-image-processing-in-2021/ for a survey of Python libraries for image processing.
Standard Python Libaries
One of the very nice aspects of image processing with Python is the wide adoption of basic NumPy arrays to represent image data. This facilitates the use of methods from multiple packages for a particular project.
NumPy NumPy provides the basic multi-dimensional array data structure used by other image processing libraries.
Matplotlib Matplotlib incudes functions to read image files in several common formats, and display images.
scipy.ndimage SciPy provides a useful collection of image processing algorithms.
scikit-image A library designed to interoperate with NumPy and SciPy.ndimage libraries.
Pure Python Libraries
Pillow A fork of the original Python image library (PIL), Pillow is one of the most important Python libraries for image processing, and the foundation for many other packages. For x86 architectures, Pillow-SIMD is a fork that "follows" Pillow to provide algorithms tuned to the x86 SIMD hardware.
Imageio Imageio is a cross-platform, pure Python library to read and write images, video, volumetric data, and scientific formats. This is a good, lightweight choice if the reading and writing files in multiple formats is the primary purpose.
Full Featured Computer Vision Packages
OpenCV/cv2 OpenCV is a C++ library launched in 1999 at Intel Research to advance real-time computer vision. In 2011 is was taken over by the non-profit foundation OpenCV.org. OpenCV is widely used with mulitple language bindings, full ecosystem of third party documentation, training courses, and books. The Python bindings are accessed by importing the cv2 Python module.
Mahotas A computer vision library written in Python and CPython. While OpenCV is the fastest package, the lack of type checking can result in hard crashes. Mahotas includes type checking with some loss of speed. A paper describing the package is available http://doi.org/10.5334/jors.ac.
SimpleITK The Natonal Library of Medicine Insight Segmentation and Registration Toolkit (ITK). book
Photo Conversion and Management Tools
ImageMagick ImageMagick is a widely-used, full-featured package to create, edit, compose, and convert images. It is available on all major platforms including the Raspberry Pi OS. Wand is a Python binding to the ImageMagick API. See https://www.pythonpool.com/imagemagick-python/ for information using Wand and ImageMagick.
XnView MP A commercial digital asset management tool available on Windows, Mac OS, and Linux. There are many alternatives on the market.
Outdated or deprecated libraries
These are included here so you know what not to use for your projects.
cImage A simple image processing library for Python. Intended for use in introductory computer science courses. Uses PIL.
PIL Python image library. Hasn't been updated since 2009, now outdated and insecure. Largely been replaced by Pillow, a friendly fork of PIL.
SimpleCV An open source Python library based on OpenCV designed for rapid development of cmputer vision applications, with accompanying book from Reilly Media. There's no evidence of continued development, and not compatible with recent versions of Python 3.
For this case study we will be developing techniques to analyze images stored as computer files. Later we may consider applying these techniques to a live video stream, but for now we'll use stored images.
Things to know about image files.
General recommendations for scientific use.
Image files also carry meta-data providing additional information about the image. There several different standards for meta data, the most common being exif data embedded in the image file, or a companion .xmp file hold data in the Extensible Metadata Platform format. exiftool is a command line tool included with many operating systems.
!/usr/local/bin/exiftool data/25-miniM.tif
ExifTool Version Number : 11.52 File Name : 25-miniM.tif Directory : data File Size : 2.4 MB File Modification Date/Time : 2021:11:23 13:25:31-05:00 File Access Date/Time : 2021:11:29 15:44:32-05:00 File Inode Change Date/Time : 2021:11:23 13:25:44-05:00 File Permissions : rw-r--r-- File Type : TIFF File Type Extension : tif MIME Type : image/tiff Exif Byte Order : Little-endian (Intel, II) Subfile Type : Full-resolution Image Image Width : 2560 Image Height : 2048 Bits Per Sample : 8 8 8 Compression : JPEG Photometric Interpretation : RGB Strip Offsets : (Binary data 1924 bytes, use -b option to extract) Orientation : Horizontal (normal) Samples Per Pixel : 3 Rows Per Strip : 8 Strip Byte Counts : (Binary data 1405 bytes, use -b option to extract) X Resolution : 300 Y Resolution : 300 Planar Configuration : Chunky Resolution Unit : inches JPEG Tables : (Binary data 289 bytes, use -b option to extract) Image Size : 2560x2048 Megapixels : 5.2
What's not present in this example? (Hint: What colors are being displayed?)
Here's the exif data for a photograph prepared through a typical photographer's workflow: RAW --> Noise Reduction --> Adobe Lightroom --> .jpg for export. Carefully look for the color space data.
!/usr/local/bin/exiftool data/eagles-photo-low-quality.jpg
ExifTool Version Number : 11.52 File Name : eagles-photo-low-quality.jpg Directory : data File Size : 248 kB File Modification Date/Time : 2021:11:29 16:25:58-05:00 File Access Date/Time : 2021:11:29 16:28:21-05:00 File Inode Change Date/Time : 2021:11:29 16:28:18-05:00 File Permissions : rw-r--r-- File Type : JPEG File Type Extension : jpg MIME Type : image/jpeg Exif Byte Order : Little-endian (Intel, II) Make : OLYMPUS CORPORATION Camera Model Name : E-M1X X Resolution : 240 Y Resolution : 240 Resolution Unit : inches Software : Adobe Lightroom 5.0 (Macintosh) Modify Date : 2021:11:29 16:25:58 Copyright : © Jeffrey Kantor Exposure Time : 1/400 F Number : 8.0 Exposure Program : Aperture-priority AE ISO : 500 Sensitivity Type : Standard Output Sensitivity Exif Version : 0231 Date/Time Original : 2021:07:15 07:54:21 Create Date : 2021:07:15 07:54:21 Offset Time : -05:00 Shutter Speed Value : 1/400 Aperture Value : 8.0 Exposure Compensation : -0.3 Max Aperture Value : 8.0 Subject Distance : 35.4 m Metering Mode : Multi-segment Light Source : Unknown Flash : Off, Did not fire Focal Length : 600.0 mm Color Space : sRGB Focal Plane X Resolution : 3010.362305 Focal Plane Y Resolution : 3010.362305 Focal Plane Resolution Unit : cm File Source : Digital Camera CFA Pattern : [Red,Green][Green,Blue] Custom Rendered : Normal Exposure Mode : Manual White Balance : Auto Digital Zoom Ratio : 1 Focal Length In 35mm Format : 1202 mm Scene Capture Type : Standard Gain Control : High gain up Contrast : Normal Saturation : Normal Sharpness : Normal Serial Number : BJ4A06646 Lens Info : 600mm f/8 Lens Model : M.300mm F4.0 + MC-20 Lens Serial Number : ACA204117 GPS Version ID : 2.2.0.0 GPS Img Direction : 134.5 Compression : JPEG (old-style) Thumbnail Offset : 1070 Thumbnail Length : 19508 Displayed Units X : inches Displayed Units Y : inches Current IPTC Digest : 786d1b99d5d40bcdba895d751f9372b7 Coded Character Set : UTF8 Application Record Version : 4 Time Created : 07:54:21 Digital Creation Date : 2021:07:15 Digital Creation Time : 07:54:21 Copyright Notice : © Jeffrey Kantor Photoshop Thumbnail : (Binary data 19508 bytes, use -b option to extract) IPTC Digest : 786d1b99d5d40bcdba895d751f9372b7 Profile CMM Type : Linotronic Profile Version : 2.1.0 Profile Class : Display Device Profile Color Space Data : RGB Profile Connection Space : XYZ Profile Date Time : 1998:02:09 06:49:00 Profile File Signature : acsp Primary Platform : Microsoft Corporation CMM Flags : Not Embedded, Independent Device Manufacturer : Hewlett-Packard Device Model : sRGB Device Attributes : Reflective, Glossy, Positive, Color Rendering Intent : Perceptual Connection Space Illuminant : 0.9642 1 0.82491 Profile Creator : Hewlett-Packard Profile ID : 0 Profile Copyright : Copyright (c) 1998 Hewlett-Packard Company Profile Description : sRGB IEC61966-2.1 Media White Point : 0.95045 1 1.08905 Media Black Point : 0 0 0 Red Matrix Column : 0.43607 0.22249 0.01392 Green Matrix Column : 0.38515 0.71687 0.09708 Blue Matrix Column : 0.14307 0.06061 0.7141 Device Mfg Desc : IEC http://www.iec.ch Device Model Desc : IEC 61966-2.1 Default RGB colour space - sRGB Viewing Cond Desc : Reference Viewing Condition in IEC61966-2.1 Viewing Cond Illuminant : 19.6445 20.3718 16.8089 Viewing Cond Surround : 3.92889 4.07439 3.36179 Viewing Cond Illuminant Type : D50 Luminance : 76.03647 80 87.12462 Measurement Observer : CIE 1931 Measurement Backing : 0 0 0 Measurement Geometry : Unknown Measurement Flare : 0.999% Measurement Illuminant : D65 Technology : Cathode Ray Tube Display Red Tone Reproduction Curve : (Binary data 2060 bytes, use -b option to extract) Green Tone Reproduction Curve : (Binary data 2060 bytes, use -b option to extract) Blue Tone Reproduction Curve : (Binary data 2060 bytes, use -b option to extract) XMP Toolkit : Adobe XMP Core 7.0-c000 1.000000, 0000/00/00-00:00:00 White Level : 29222 Adobe White Level : 29944 Creator Tool : Adobe Lightroom 5.0 (Macintosh) Rating : 4 Metadata Date : 2021:11:29 16:25:58-05:00 Lens : M.300mm F4.0 + MC-20 Approximate Focus Distance : 35.4 Flash Compensation : 0 Date Created : 2021:07:15 07:54:21 Document ID : xmp.did:35e88d02-dcad-40ee-954c-d74068c1e965 Original Document ID : DD64D26CAFF63683666FACABB6CE7C1D Instance ID : xmp.iid:35e88d02-dcad-40ee-954c-d74068c1e965 Format : image/jpeg Raw File Name : Version : 14.0 Process Version : 11.0 Color Temperature : 4926 Tint : +15 Exposure 2012 : +0.03 Contrast 2012 : +33 Highlights 2012 : -66 Shadows 2012 : +5 Whites 2012 : +2 Blacks 2012 : -21 Texture : +40 Clarity 2012 : +13 Dehaze : +8 Vibrance : +20 Parametric Shadows : 0 Parametric Darks : 0 Parametric Lights : 0 Parametric Highlights : 0 Parametric Shadow Split : 25 Parametric Midtone Split : 50 Parametric Highlight Split : 75 Sharpen Radius : +1.0 Sharpen Detail : 25 Sharpen Edge Masking : 0 Luminance Smoothing : 0 Color Noise Reduction : 0 Hue Adjustment Red : 0 Hue Adjustment Orange : 0 Hue Adjustment Yellow : 0 Hue Adjustment Green : 0 Hue Adjustment Aqua : 0 Hue Adjustment Blue : 0 Hue Adjustment Purple : 0 Hue Adjustment Magenta : 0 Saturation Adjustment Red : 0 Saturation Adjustment Orange : 0 Saturation Adjustment Yellow : 0 Saturation Adjustment Green : 0 Saturation Adjustment Aqua : 0 Saturation Adjustment Blue : 0 Saturation Adjustment Purple : 0 Saturation Adjustment Magenta : 0 Luminance Adjustment Red : 0 Luminance Adjustment Orange : 0 Luminance Adjustment Yellow : 0 Luminance Adjustment Green : 0 Luminance Adjustment Aqua : 0 Luminance Adjustment Blue : 0 Luminance Adjustment Purple : 0 Luminance Adjustment Magenta : 0 Split Toning Shadow Hue : 0 Split Toning Shadow Saturation : 0 Split Toning Highlight Hue : 0 Split Toning Highlight Saturation: 0 Split Toning Balance : 0 Color Grade Midtone Hue : 0 Color Grade Midtone Sat : 0 Color Grade Shadow Lum : 0 Color Grade Midtone Lum : 0 Color Grade Highlight Lum : 0 Color Grade Blending : 50 Color Grade Global Hue : 0 Color Grade Global Sat : 0 Color Grade Global Lum : 0 Auto Lateral CA : 0 Lens Profile Enable : 0 Lens Manual Distortion Amount : 0 Vignette Amount : 0 Defringe Purple Amount : 0 Defringe Purple Hue Lo : 30 Defringe Purple Hue Hi : 70 Defringe Green Amount : 0 Defringe Green Hue Lo : 40 Defringe Green Hue Hi : 60 Perspective Upright : 0 Perspective Vertical : 0 Perspective Horizontal : 0 Perspective Rotate : 0.0 Perspective Aspect : 0 Perspective Scale : 100 Perspective X : 0.00 Perspective Y : 0.00 Grain Amount : 0 Post Crop Vignette Amount : 0 Shadow Tint : 0 Red Hue : 0 Red Saturation : 0 Green Hue : 0 Green Saturation : 0 Blue Hue : 0 Blue Saturation : 0 Override Look Vignette : False Tone Curve Name 2012 : Linear Camera Profile : Adobe Standard Camera Profile Digest : 17201179F5294A1E72C1552BAE550EEF Auto Tone Digest : 1715A4BED36C368980DA388E7373CCEE Auto Tone Digest No Sat : 3E753A3B412B959232482F410AEDEFE8 Has Settings : True Crop Top : 0.078092 Crop Left : 0.170981 Crop Bottom : 0.956094 Crop Right : 0.957141 Crop Angle : 0.192816 Crop Constrain To Warp : 0 Has Crop : True Already Applied : True History Action : edited, derived, saved History When : 2021:07:15 22:11:47Z, 2021:11:29 16:25:58-05:00 History Parameters : converted from image/dng to image/jpeg, saved to new location History Instance ID : xmp.iid:35e88d02-dcad-40ee-954c-d74068c1e965 History Software Agent : Adobe Lightroom 5.0 (Macintosh) History Changed : / Derived From Document ID : DD64D26CAFF63683666FACABB6CE7C1D Derived From Original Document ID: DD64D26CAFF63683666FACABB6CE7C1D Rights : © Jeffrey Kantor Tone Curve PV2012 : 0, 0, 255, 255 Tone Curve PV2012 Red : 0, 0, 255, 255 Tone Curve PV2012 Green : 0, 0, 255, 255 Tone Curve PV2012 Blue : 0, 0, 255, 255 Look Name : Adobe Portrait Look Amount : 1 Look Uuid : D6496412E06A83789C499DF9540AA616 Look Supports Amount : false Look Supports Monochrome : false Look Supports Output Referred : false Look Group : Profiles Look Parameters Version : 14.0 Look Parameters Process Version : 11.0 Look Parameters Convert To Grayscale: False Look Parameters Camera Profile : Adobe Standard Look Parameters Look Table : E5A76DBB8B3F132A04C01AF45DC2EF1B Look Parameters Tone Curve PV2012: 0, 0, 66, 64, 190, 192, 255, 255 DCT Encode Version : 100 APP14 Flags 0 : [14], Encoded with Blend=1 downsampling APP14 Flags 1 : (none) Color Transform : YCbCr Image Width : 1704 Image Height : 2048 Encoding Process : Baseline DCT, Huffman coding Bits Per Sample : 8 Color Components : 3 Y Cb Cr Sub Sampling : YCbCr4:2:0 (2 2) Aperture : 8.0 Date/Time Created : 2021:07:15 07:54:21 Digital Creation Date/Time : 2021:07:15 07:54:21 Image Size : 1704x2048 Megapixels : 3.5 Scale Factor To 35 mm Equivalent: 2.0 Shutter Speed : 1/400 Modify Date : 2021:11:29 16:25:58-05:00 Thumbnail Image : (Binary data 19508 bytes, use -b option to extract) Circle Of Confusion : 0.015 mm Depth Of Field : 0.82 m (34.99 - 35.82 m) Field Of View : 1.7 deg Focal Length : 600.0 mm (35 mm equivalent: 1202.0 mm) Hyperfocal Distance : 3000.37 m Light Value : 12.3
Accurate color reproduction requires color management from source to display. At each Color management is now a standard part of most operating systems
Our approach ...
We track overall code dependencies by consolidating imports into this cell. Note that we'll be using elements from multiple packages by relying on the underlying NumPy representation of images to hold the current state of the process.
%matplotlib inline
# standard Python libraries
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
from scipy import ndimage
# computer vision libraries
import cv2 as cv2
import mahotas as mh
As a first step, read the image, convert to rgb scale, and display. All of these packages have a means of reading raster file images in common formats. There are small (and sometimes frustrating) differences among them. Here we use the Matplotlib imread() method which reads and returns a numpy array.
The array will typical have 2 or 3 dimensions $(h, w, d)$ where $h$ and $w$ are image height and width, and $d$ is pixel depth.
filepath = "data/25-miniM.tif"
# read color image with OpenCV
img_bgr = cv2.imread(filepath)
print(img_bgr.shape)
# convert to RGB
img_rgb = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2RGB)
# display images with Matploblib
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].imshow(img_rgb)
ax[0].set_title("RGB image displayed as RGB")
ax[1].imshow(img_bgr)
ax[1].set_title("A BGR image incorrectly shown as RGB")
(2048, 2560, 3)
Text(0.5, 1.0, 'A BGR image incorrectly shown as RGB')
# reading files with Pillow
img = Image.open(filepath)
plt.imshow(img)
<matplotlib.image.AxesImage at 0x7feb9a93e5e0>
Observations:
crop_img = img_rgb[100:2000, 100:2400, :].copy()
plt.imshow(crop_img)
<matplotlib.image.AxesImage at 0x7feb04e3b490>
An image is comprised of one or more channels
r, g, b = cv2.split(crop_img)
Histograms are a tool for analyzing the distribution of gray levels in a channel. It's a powerful tool for controlling exposure and processing images for presentation.
def histogram(channel, bp=3, wp=252):
"""Return histogram and bins for a single channel."""
hist = cv2.calcHist([channel], [0], None, [wp-bp+1], [bp, wp])
bins = np.array([b for b in range(bp, wp+1)])
return hist.flatten(), bins
def display_channel(channel, label=""):
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ax[0].imshow(channel.astype(np.uint8), cmap="gray")
ax[0].set_title(label)
hist, bins = histogram(channel.astype(np.uint8))
ax[1].fill_between(bins, hist, alpha=0.4, color="k", label=label)
ax[1].legend()
display_channel(r, "red")
display_channel(g, "green")
display_channel(b, "blue")
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ax[0].imshow(img)
for color, channel in zip(['r', 'g', 'b'], [r, g, b]):
hist, bins = histogram(channel)
ax[1].fill_between(bins, hist, color=color, label=color, alpha=0.3)
ax[1].legend()
<matplotlib.legend.Legend at 0x7feb05e495e0>
We see the blue leds used to excite the flourophores bleed over into the green channel. It would be best if this could be corrected in the experiment, perhaps by positioning a bandpass filter in front of the leds. What we will attempt here is subtract a multiple of blue channel from the green channel, followed by exposure adjustments. The goal is to provide a cleaner image for doing particle labeling and counting.
By trial and error, we find a weighted difference of the green and blue channels, and a rescaling of the tone curve that retains the particles and reduces background interference.
blue_weight = 0.58
offset = -20
# subtract blue channel from green channel
cimg = (g - blue_weight*b) + offset
display_channel(cimg, "particle channel'")
# set zero threshold and convert to integer
cimg = np.where(cimg < 0, 0, cimg)
cimg = cimg.astype(np.uint8)
display_channel(cimg, "particle channel zero threshold")
himg = cv2.equalizeHist(cimg)
display_channel(himg, 'equalized')
Observations
# kernel size (odd number)
ksize = 21
# median filter
bimg = cv2.medianBlur(himg, ksize)
display_channel(bimg, "Low Pass Filter")
# Gaussian filter
bimg = cv2.GaussianBlur(himg, (ksize, ksize), 0)
display_channel(bimg, "Low Pass Filter")
The purpose of threshold is to isolate the features of interest from background noise.
https://docs.opencv.org/3.4/d7/d4d/tutorial_py_thresholding.html
threshold = 130
T, simg = cv2.threshold(bimg, threshold, 255, cv2.THRESH_BINARY)
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ax[0].imshow(crop_img)
ax[1].imshow(simg, cmap="gray")
<matplotlib.image.AxesImage at 0x7feb0c643580>
The next goal is to remove noise and to separate particles.
kernel = np.ones((3, 3))
iterations = 8
dimg = cv2.dilate(simg, kernel=kernel, iterations=iterations)
eimg = cv2.erode(dimg, kernel=kernel, iterations=iterations)
fig, ax =
plt.imshow(eimg, cmap="gray")
iterations = 18
fimg = cv2.erode(eimg, kernel=kernel, iterations=iterations)
gimg = cv2.dilate(fimg, kernel=kernel, iterations=iterations)
plt.imshow(gimg, cmap="gray")
<matplotlib.image.AxesImage at 0x7feb13bb32e0>
mimg = cv2.morphologyEx(simg, cv2.MORPH_CLOSE, kernel=kernel, iterations=7)
ming = cv2.morphologyEx(mimg, cv2.MORPH_OPEN, kernel=kernel, iterations=50)
plt.imshow(mimg, cmap="gray")
<matplotlib.image.AxesImage at 0x7feb044bf3a0>
http://pageperso.lif.univ-mrs.fr/~francois.denis/IAAM1/scipy-html-1.0.0/tutorial/ndimage.html
structure = np.ones((3, 3))
# find particles and plot objects
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
label_img, particle_count = ndimage.measurements.label(mimg, structure=structure)
# pixels are labeled with numbers corresponding to objects
ax.imshow(np.where(label_img > 0, 1, np.zeros(mimg.shape)), cmap="gray")
# centers of mass
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
pts = ndimage.center_of_mass(mimg, label_img, np.arange(1, particle_count + 1))
ax[0].imshow(img)
x = [x for (y,x) in pts]
y = [y for (y,x) in pts]
ax[0].plot(x, y, '.', ms=5, color='r')
[<matplotlib.lines.Line2D at 0x7feb02cdd0a0>]
# find the objects and place bounding boxes
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
slices = ndimage.find_objects(label_img)
sizes = [np.sqrt((b.stop-b.start)**2 + (a.stop-a.start)**2) for a,b in slices]
ax[1].hist(sizes, bins=200)
ax[1].set_xlim(0, 200)
ax[0].imshow(img)
for k,size in enumerate(sizes):
if size > 50 and size < 200:
a,b = slices[k]
ya = a.start
yb = a.stop
xa = b.start
xb = b.stop
ax[0].plot([xa, xb, xb, xa, xa], [yb, yb, ya, ya, yb], 'w')
We need to improve the quality of the image
Image processing steps
To facilitate embedded use in a device, the next step is to consolidate these procedures into a class.
import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv
from scipy import ndimage
class Channel():
def __init__(self, data):
self.data = data
@property
def histogram(self):
return cv.calcHist([self.img], [0], None, [256], [0, 255])
def imshow(self, ax=None, title=None):
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.imshow(self.img, cmap="gray")
if not title is None:
ax.set_title(title)
def display_histogram(self, ax, color='k'):
bins = np.arange(0, 256)
ax.plot(bins, self.histogram, color=color)
class Particle():
def __init__(self, x, y, slice):
self.x = x
self.y = y
self.slice = slice
a, b = slice
self.xmin = b.start
self.ymin = a.start
self.xmax = b.stop
self.ymax = a.stop
def size(self):
return np.sqrt((self.xmax - self.xmin) * (self.ymax - self.ymin))
def bounding_box(self):
return (self.xmin, self.xmax, self.xmax, self.xmin, self.xmin), \
(self.ymin, self.ymin, self.ymax, self.ymax, self.ymin)
class ParticleLabeler():
def __init__(self, filepath=None):
self.filepath = filepath
self.channels = dict()
self.particles = list()
if not filepath is None:
self.imread(filepath)
@property
def particle_count(self):
return len(self.particles)
def imread(self, filepath):
"""Read image from given filepath."""
self.filepath = filepath
b, g, r = cv.split(cv.imread(filepath))
self.channels['r'] = Channel(r)
self.channels['g'] = Channel(g)
self.channels['b'] = Channel(b)
def imshow(self, ax=None, rgb=('r','g','b')):
"""Display image on given plot axis using specified rgb channels."""
r, g, b = rgb
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.imshow(np.dstack((self.channels[r].img,
self.channels[g].img,
self.channels[b].img)))
def display_histogram(self, ax=None):
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
for color, channel in self.channels.items():
channel.display_histogram(ax, color=color)
def composite(self, rgb_weights={'r':1, 'g':1, 'b':1}, dst='k'):
composite = sum([rgb_weights[c]*self.channels[c].img for c in rgb_weights.keys()])
composite = np.where(composite < 0, 0, composite)
self.channels[dst] = Channel(composite.astype(np.uint8))
def threshold(self, src='k', dst='k', r=71):
blurred = cv.medianBlur(self.channels[src].img, r)
T, img = cv.threshold(blurred, 0, 255, cv.THRESH_OTSU)
self.channels[dst] = Channel(img)
def adaptive_threshold(self, src='k', dst='k'):
blurred = cv.medianBlur(self.channels[src].img, 5)
img = cv.adaptiveThreshold(blurred, 255, cv.ADAPTIVE_THRESH_GAUSSIAN_C, cv.THRESH_BINARY, 101, 10)
self.channels[dst] = Channel(img)
def label(self, src='k'):
img = self.channels[src].img
labels, particle_count = ndimage.measurements.label(img, structure=np.ones((3, 3)))
yx_pts = ndimage.center_of_mass(img, labels, np.arange(1, particle_count + 1))
slices = ndimage.find_objects(labels)
self.particles = [Particle(yx[1], yx[0], slice) for yx, slice in zip(yx_pts, slices)]
def find_particles(self, filepath, size=(50, 150), ax=None):
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(10, 12))
self.imread(filepath)
self.composite({'g':1, 'b':-0.53}, 'k')
self.threshold(src='k', dst='k')
self.label(src='k')
self.imshow(ax)
k = 0
for p in labeler.particles:
if p.size() > 50 and p.size() < 150:
k += 1
x, y = p.bounding_box()
ax.plot(x, y, 'w')
ax.plot([p.x], [p.y], 'r.', ms=5)
ax.text(p.x + 40, p.y - 40, f"{k}", color='white')
help(Channel)
Help on class Channel in module __main__: class Channel(builtins.object) | Channel(img) | | A class to represent one channel of a multi-channel image. | ... | | Attributes | ---------- | img : | | Methods defined here: | | __init__(self, img) | Initialize self. See help(type(self)) for accurate signature. | | display_histogram(self, ax, color='k') | | imshow(self, ax=None, title=None) | | ---------------------------------------------------------------------- | Readonly properties defined here: | | histogram | | ---------------------------------------------------------------------- | Data descriptors defined here: | | __dict__ | dictionary for instance variables (if defined) | | __weakref__ | list of weak references to the object (if defined)
# TESTING
# create labeler instance
labeler = ParticleLabeler()
# read image file
labeler.imread("data/25-miniM.tif")
# display image and histogram
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
labeler.imshow(ax[0])
labeler.display_histogram(ax[1])
ax[1].set_ylim(0, 150000)
# make composite channel
labeler.composite({'g':1, 'b':-0.53}, 'k')
# threshold
labeler.threshold(src='k', dst='k')
# label
labeler.label(src='k')
print(labeler.particle_count)
# display all channels
for color, channel in labeler.channels.items():
fig, ax = plt.subplots(1, 2, figsize=(15,5))
channel.imshow(ax[0], color)
channel.display_histogram(ax[1], color=color)
ax[1].set_ylim(0, 150000)
fix, ax = plt.subplots(1, 1, figsize=(12,10))
labeler.imshow(ax)
for p in labeler.particles:
if p.size() > 50 and p.size() < 150:
x, y = p.bounding_box()
ax.plot(x, y, 'w')
ax.plot([p.x], [p.y], 'r.', ms=5)
46
labeler = ParticleLabeler()
labeler.find_particles("data/25-miniM.tif", size=(50, 150))
print(labeler.particle_count)
46
img = labeler.channels['g'].img
fimg = cv.GaussianBlur(img, (71, 71), 0)
fimg = cv.medianBlur(img, 51)
T, timg = cv.threshold(fimg, 80, 255, cv.THRESH_BINARY)
#T, timg = cv.threshold(fimg, 0, 255, cv.THRESH_OTSU)
fig, ax = plt.subplots(1, 1, figsize=(10, 12))
#ax.imshow(img, cmap="gray")
#ax.imshow(img, cmap="gray")
ax.imshow(timg, cmap="gray")
ax.set_title(f"Threshold = {T}")
#ax.imshow(edges, cmap="hot", alpha=0.8)
Text(0.5, 1.0, 'Threshold = 80.0')
# TESTING
# create labeler instance
labeler = ParticleLabeler()
# read image file
labeler.imread("data/25-camera.jpg")
# make composite channel
labeler.composite({'g':1}, 'k')
# threshold
labeler.threshold(src='k', dst='k', r=41)
# label
labeler.label(src='k')
print(labeler.particle_count)
fix, ax = plt.subplots(3, 1, figsize=(12,20))
labeler.imshow(ax[0])
labeler.channels['k'].imshow(ax[1])
labeler.imshow(ax[2])
for p in labeler.particles:
if p.size() > 50 and p.size() < 150:
x, y = p.bounding_box()
ax[2].plot(x, y, 'w')
ax[2].plot([p.x], [p.y], 'r.', ms=5)
34
labeler = ParticleLabeler()
labeler.imread("data/25-camera.jpg")
labeler.imshow()
labeler.channels['g'].imshow(title="Green Channel")
labeler.composite({'g':1, 'b':-0.53}, 'k')
labeler.channels['k'].imshow(title="Composite Green - Blue")
labeler.adaptive_threshold(src='g', dst='k2')
labeler.channels['k2'].imshow()
circles = cv.HoughCircles(timg,
cv.HOUGH_GRADIENT,
1.5,
30,
param1=80,
param2=20,
minRadius=50,
maxRadius=200)
fig, ax = plt.subplots(1, 1, figsize=(10, 12))
ax.imshow(img, cmap="gray")
for circle in circles[0]:
x, y, r = circle
plt.plot([x], [y], 'r.', ms=20)
